clear all

M = dlmread('usawk.csv'); %loading KLEMS data

%1: year, 1947-2014
%2: industry, 1-65
%3: current Y
%4: current K
%5: current L
%6: current interm
%7: constant Y
%8: constant K
%9: constant L
%10: constant interm

T=2014-1947+1;
S=61;
vecS=[1:61]; %selecting industries, excluding Federal General government, Federal Government enterprises, S&L Government enterprises, S&L General Government

for t=1:T
    for s=vecS;
        
        matY_curr(t,s)=M((s-1)*T+t,3);
        matK_curr(t,s)=M((s-1)*T+t,4);
        matL_curr(t,s)=M((s-1)*T+t,5);
        matI_curr(t,s)=M((s-1)*T+t,6);
        matVA_curr(t,s)=matY_curr(t,s)-matI_curr(t,s);
                
        matY_const(t,s)=M((s-1)*T+t,7);
        matK_const(t,s)=M((s-1)*T+t,8);
        matL_const(t,s)=M((s-1)*T+t,9);
        matI_const(t,s)=M((s-1)*T+t,10);
               
        matsK(t,s)=matK_curr(t,s)/matY_curr(t,s);
        matsL(t,s)=matL_curr(t,s)/matY_curr(t,s);
        matsI(t,s)=matI_curr(t,s)/matY_curr(t,s);
        
    end
    
    vecGDP_curr(t,1)=sum(matVA_curr(t,:));
    vecK_curr(t,1)=sum(matK_curr(t,:));
    vecL_curr(t,1)=sum(matL_curr(t,:));
    
    %current domar weights
    for s=vecS;
        matlambda(t,s)=matY_curr(t,s)/vecGDP_curr(t); 
        matlambdaK(t,s)=matK_curr(t,s)/vecK_curr(t); 
        matlambdaL(t,s)=matL_curr(t,s)/vecL_curr(t); 
    end
end


for t=1:T-1
    for s=vecS;        
        %average domar weights and average input shares 
        matavlambda(t,s)=0.5*(matlambda(t,s)+ matlambda(t+1,s)); 
        matavlambdaK(t,s)=0.5*(matlambdaK(t,s)+ matlambdaK(t+1,s)); 
        matavlambdaL(t,s)=0.5*(matlambdaL(t,s)+ matlambdaL(t+1,s)); 
        matavsK(t,s)=0.5*(matsK(t,s)+matsK(t+1,s));
        matavsL(t,s)=0.5*(matsL(t,s)+matsL(t+1,s));
        matavsI(t,s)=0.5*(matsI(t,s)+matsI(t+1,s));
            
        %growth rates
        matgY(t,s)=log(matY_const(t+1,s))-log(matY_const(t,s));
        matgK(t,s)=log(matK_const(t+1,s))-log(matK_const(t,s));
        matgL(t,s)=log(matL_const(t+1,s))-log(matL_const(t,s));
        matgI(t,s)=log(matI_const(t+1,s))-log(matI_const(t,s));
            
        %using averages
        matgTFPav(t,s)=matgY(t,s)-matavsK(t,s)*matgK(t,s)-matavsL(t,s)*matgL(t,s)-matavsI(t,s)*matgI(t,s);
        
    end
end


for t=1:T-1
        
    vecavTFP_actual(t,1)=sum(matgTFPav(t,:).*matavlambda(t,:)); %only substitution effects
    vecavTFP_init(t,1)=sum(matgTFPav(t,:).*matavlambda(1,:)); %intial shares
    TFP_EV(t,1)=sum(sum(matgTFPav(1:t,:),1).*matavlambda(t,:)); %only taste / income effects

end


TFP_avactual=cumsum(vecavTFP_actual); %only substitution effects
TFP_avinit=cumsum(vecavTFP_init); %intial shares

plot([1948:2014],([TFP_avactual TFP_EV TFP_avinit]),'LineWidth',2); legend('Welfare with only substitution effect','Welfare with only taste/income effect','Initial Shares','Location','NorthWest')